Econometría II: ARIMA

Departamento de Economía

Carlos A. Yanes G.

2023-09-06

Paquetes con que se trabaja la sesión

R libraries

Los paquetes que vamos a utilizar y que se deben instalar para ser usados en clase son:

install.packages(c("pacman","TSstudio", "urca", "forecast", "devtools"))

# Si tiene problemas con TSstudio podemos mirar entonces en Github
devtools::install_github("RamiKrispin/TSstudio")

Ejecución

No olvide que debe desde luego tener presente en su documento de R Markdown cargar estos paquetes.

library(pacman)
p_load(tidyverse, TSstudio, tseries, urca, forecast)

Preambulo

Modelar ARIMA

Debemos primero mirar la estacionariedad.
Si no es estacionaria la serie hay que diferenciar.
Después de testear procedemos a leer sus funciones de autocovarianzas denominadas correlogramas
Constatar el RUIDO BLANCO
Predecir

Datos

Cartera comercial de bancos

Enero de 2016 hasta Abril 2023
Web site de Superintendencia Financiera
Frecuencia Mensual

Agradecimientos

Esta sesión se construye a partir de los elementos de los paquetes de tseries y (Krispin 2019), quien desarrolló el paquete TSstudio respectivamente

A lo que vinimos

Cargar datos

library(pacman)
p_load(readxl, TSstudio, tidyverse, stats, urca, forecast, ggfortify, ggplot2, tseries)
bd <- read_excel("cartera.xlsx")
bd <- bd |> select(Cartera)
cartera <- ts(bd, frequency=12, start=c(2016,1))

Graficar normal

Code

ts_plot(cartera,
        title = "Cartera comercial de Bancos", 
        Ytitle = "En Miles de Millones de $",
        Xtitle = "Años")

Code

prim <- ur.df(y=cartera, lags=0, type='none')
summary(prim)


############################################### 
# Augmented Dickey-Fuller Test Unit Root Test # 
############################################### 

Test regression none 


Call:
lm(formula = z.diff ~ z.lag.1 - 1)

Residuals:
   Min     1Q Median     3Q    Max 
 -8614  -1856   -262   1731  14834 

Coefficients:
         Estimate Std. Error t value Pr(>|t|)    
z.lag.1 0.0068402  0.0008212   8.329 1.14e-12 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3653 on 86 degrees of freedom
Multiple R-squared:  0.4465,    Adjusted R-squared:  0.4401 
F-statistic: 69.38 on 1 and 86 DF,  p-value: 1.139e-12


Value of test-statistic is: 8.3292 

Critical values for test statistics: 
     1pct  5pct 10pct
tau1 -2.6 -1.95 -1.61

Code

prim2 <- ur.df(y=cartera, lags=0, type='drift')
summary(prim2)


############################################### 
# Augmented Dickey-Fuller Test Unit Root Test # 
############################################### 

Test regression drift 


Call:
lm(formula = z.diff ~ z.lag.1 + 1)

Residuals:
    Min      1Q  Median      3Q     Max 
-8991.9 -1598.5   365.7  1936.5 14876.5 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)   
(Intercept) -4.283e+03  2.453e+03  -1.746  0.08443 . 
z.lag.1      1.571e-02  5.144e-03   3.054  0.00302 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3610 on 85 degrees of freedom
Multiple R-squared:  0.09887,   Adjusted R-squared:  0.08827 
F-statistic: 9.326 on 1 and 85 DF,  p-value: 0.003017


Value of test-statistic is: 3.0538 37.0386 

Critical values for test statistics: 
      1pct  5pct 10pct
tau2 -3.51 -2.89 -2.58
phi1  6.70  4.71  3.86

Code

prim3 <- ur.df(y=cartera, lags=0, type='trend')
summary(prim3)


############################################### 
# Augmented Dickey-Fuller Test Unit Root Test # 
############################################### 

Test regression trend 


Call:
lm(formula = z.diff ~ z.lag.1 + 1 + tt)

Residuals:
    Min      1Q  Median      3Q     Max 
-8860.2 -1405.5   318.1  1848.6 14698.6 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept) 8.965e+02  6.550e+03   0.137    0.891
z.lag.1     1.964e-04  1.890e-02   0.010    0.992
tt          4.831e+01  5.663e+01   0.853    0.396

Residual standard error: 3616 on 84 degrees of freedom
Multiple R-squared:  0.1066,    Adjusted R-squared:  0.08534 
F-statistic: 5.012 on 2 and 84 DF,  p-value: 0.008785


Value of test-statistic is: 0.0104 24.856 5.0119 

Critical values for test statistics: 
      1pct  5pct 10pct
tau3 -4.04 -3.45 -3.15
phi2  6.50  4.88  4.16
phi3  8.73  6.49  5.47

Code

# adf.test(cartera,k=0)
adf.test(cartera) # automático


    Augmented Dickey-Fuller Test

data:  cartera
Dickey-Fuller = -1.6971, Lag order = 4, p-value = 0.7006
alternative hypothesis: stationary

Analisis de funciones de autocorrelación

Anotaciones
ACF
PAC

Permiten mirar o -sospechar- si la serie es estacionaria
La serie de cartera por ninguna parte parece ser estacionaria

Code

acf(cartera, lag.max = 20)

Code

pacf(cartera, lag.max = 20)

Comportamiento Lags

Code

ts_lags(cartera)

Comportamiento Lags

Code

ts_lags(cartera, lags=c(1,4,6,11))

Diferenciemos

Code

ts_plot(diff(cartera, lag = 1),
        title = "Transformación en Nivel",
        Xtitle = "Años",
        Ytitle = "Diferencia de cartera")

Diferenciación opcional

Code

ts_plot(diff(log(cartera), lag = 1),
        title = "Transformación en logaritmo",
        Xtitle = "Años",
        Ytitle = "Diferencia/logaritmica de cartera")

Nuevamente test e identificación

DF+
DF C+T
Recuerde

Code

dcartera<-diff(cartera, lag=1) # creo el objeto
adf.test(dcartera)


    Augmented Dickey-Fuller Test

data:  dcartera
Dickey-Fuller = -3.0574, Lag order = 4, p-value = 0.1415
alternative hypothesis: stationary

Code

prim4 <- ur.df(dcartera, lags=0, type='trend')
summary(prim4)


############################################### 
# Augmented Dickey-Fuller Test Unit Root Test # 
############################################### 

Test regression trend 


Call:
lm(formula = z.diff ~ z.lag.1 + 1 + tt)

Residuals:
    Min      1Q  Median      3Q     Max 
-8738.8 -1752.5   150.9  1714.3 13670.4 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 512.7969   731.5826   0.701   0.4853    
z.lag.1      -0.6126     0.1010  -6.066 3.74e-08 ***
tt           31.7379    15.3805   2.064   0.0422 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3336 on 83 degrees of freedom
Multiple R-squared:  0.3072,    Adjusted R-squared:  0.2905 
F-statistic:  18.4 on 2 and 83 DF,  p-value: 2.433e-07


Value of test-statistic is: -6.0661 12.2669 18.3991 

Critical values for test statistics: 
      1pct  5pct 10pct
tau3 -4.04 -3.45 -3.15
phi2  6.50  4.88  4.16
phi3  8.73  6.49  5.47

Identificación para pronostico

Code

par(mfrow=c(1,2)) 
acf(dcartera)
pacf(dcartera)

Modelos

Modelo 1
Modelo 2

Code

arima_1 <- arima(cartera, order = c(1,1,0))
arima_1


Call:
arima(x = cartera, order = c(1, 1, 0))

Coefficients:
         ar1
      0.6728
s.e.  0.0782

sigma^2 estimated as 12927179:  log likelihood = -836.05,  aic = 1676.11

Code

arima_2 <- arima(cartera, order = c(2,1,0))
arima_2


Call:
arima(x = cartera, order = c(2, 1, 0))

Coefficients:
         ar1     ar2
      0.4983  0.2529
s.e.  0.1041  0.1040

sigma^2 estimated as 12090598:  log likelihood = -833.2,  aic = 1672.41

Tenemos modelo!!

Test de Ruido Blanco

Code

forecast::checkresiduals(arima(cartera, order = c(1, 1, 0), include.mean = FALSE))


    Ljung-Box test

data:  Residuals from ARIMA(1,1,0)
Q* = 37.168, df = 17, p-value = 0.003193

Model df: 1.   Total lags used: 18

Pronostico de cartera

Code

arima_fc <- forecast(arima_1, h = 12)
plot_forecast(arima_fc,
 title = "Pronostico modelo AR(1)", 
 Ytitle = "Valores",
 Xtitle = "Años")

Detras de la serie

Code

# Crear el gráfico de las dos series de tiempo
plot.ts(cartera, col = "blue", xlab = "Tiempo", ylab = "Cartera Comercial", main = "Predicción de cartera comercial")
lines(arima_fc$fitted, col = "red")
legend("topleft", legend = c("Serie Original", "Predicción"), col = c("blue", "red"), lty = 1)

Valores

Code

arima_fc

         Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
May 2023       640573.6 635965.9 645181.4 633526.7 647620.6
Jun 2023       641772.5 632792.5 650752.5 628038.7 655506.2
Jul 2023       642579.0 629291.8 655866.3 622258.0 662900.1
Aug 2023       643121.7 625746.0 660497.3 616547.9 669695.4
Sep 2023       643486.8 622289.9 664683.6 611069.0 675904.5
Oct 2023       643732.4 618983.1 668481.6 605881.7 681583.1
Nov 2023       643897.6 615846.9 671948.4 600997.7 686797.5
Dec 2023       644008.8 612882.6 675134.9 596405.5 691612.1
Jan 2024       644083.6 610082.3 678084.8 592083.1 696084.0
Feb 2024       644133.9 607433.6 680834.2 588005.6 700262.2
Mar 2024       644167.8 604922.8 683412.7 584147.8 704187.7
Apr 2024       644190.5 602536.5 685844.6 580486.2 707894.9

Selección de modelos

Existen varias formas de buscar elegir un modelo entre varios que se han estimado por pensar que es nuestro mejor modelo generador de datos. Estos vienen a ser los criterios de Akaike (AIC) y Bayes (BIC)

\(AIC=-2\times log(L)+ 2\times p\)
P es el número de parámetros del modelo
L la razón de máxima verosimilitud del modelo

\(BIC=-2\times log(L)+ p \times log(n)\)
n es es tamaño de la muestra

Gracias por su atención!!

Referencias

Krispin, Rami. 2019. Hands-on Time Series Analysis with r: Perform Time Series Analysis and Forecasting Using r. Packt Publishing Ltd.